#include<poisson_construct.h>



double fx_dot_fx[5]={1/120.0,13/60.0,11/20.0,13/60.0,1/120.0},fx_dot_Dfx[5]={-1/24.0,-5/12.0,0, 5/12.0,1/24.0}, fx_dot_DDfx[5]={1/6.0,1/3.0, -1,1/3.0,1/6.0};
//fx_dot_DDfxi=DDfx_dot_fxi
//fx_dot_Dfxi=-Dfx_dot_fxi

void poisson_grid_free(Poisson_Grid*pg)
{
    if(pg==NULL)
    {
        return;
    }
    if(pg->pgps!=NULL&&pg->size>0)
    {
        for(int i=0;i<pg->size;i++)
        {
            for(int j=0;j<pg->size;j++)
            {
                free(pg->pgps[i][j]); 

            }
            free(pg->pgps[i]);
        }
        free(pg->pgps);
    }
    free(pg);
}
void poisson_grid_init_with_depth(Poisson_Grid*pg,int depth)
{
    pg->depth=depth;
    pg->size=pow(2,depth);
    pg->width=2*pg->range/(double)(pg->size);
    double *pos_xs=(double*)malloc(sizeof(double)*pg->size);
    memset(pos_xs,0,sizeof(double)*pg->size);
    for(int i=0;i<pg->size;i++)
    {
        pos_xs[i]=(i+0.5)*pg->width-pg->range;
    }
    pg->pgps=(Poisson_Grid_Prop***)malloc(sizeof(Poisson_Grid_Prop**)*pg->size); 
    for(int i=0;i<pg->size;i++)
    {
        pg->pgps[i]=(Poisson_Grid_Prop**)malloc(sizeof(Poisson_Grid_Prop*)*pg->size);
        for(int j=0;j<pg->size;j++)
        {
            pg->pgps[i][j]=(Poisson_Grid_Prop*)malloc(sizeof(Poisson_Grid_Prop)*pg->size);
            for(int k=0;k<pg->size;k++)
            {
                poisson_grid_prop_init(&(pg->pgps[i][j][k])); 
                pg->pgps[i][j][k].pos[0]=pos_xs[i];
                pg->pgps[i][j][k].pos[1]=pos_xs[j];
                pg->pgps[i][j][k].pos[2]=pos_xs[k];
            }
        }
    }
    free(pos_xs);
}



double* poisson_compute_b(Poisson_Grid*pg )
{
    double *re=(double*)malloc(sizeof(double)*pg->size*pg->size*pg->size);
    memset(re,0,sizeof(double)*pg->size*pg->size*pg->size);
    int sum=0,id=0;
    for(int i=0;i<pg->size;i++)
    {
        for(int j=0;j<pg->size;j++)
        {
            for(int k=0;k<pg->size;k++)
            {
                id=threeD_index2_id(i,j,k,pg->size);
                for(int l=0;l<5;l++ )
                {
                    for(int m=0;m<5;m++)
                    {
                        for(int n=0;n<5;n++)
                        {
                            if(i+l-2<0||i+l-2>=pg->size||j+m-2<0||j+m-2>=pg->size||k+n-2<0||k+n-2>=pg->size)
                            {
                                continue;
                            }
                            re[id]+= pg->pgps[i+l-2][j+m-2][k+n-2].normal_coeff[0]*fx_dot_Dfx[l]*fx_dot_fx[m]*fx_dot_fx[n];
                            re[id]+= pg->pgps[i+l-2][j+m-2][k+n-2].normal_coeff[1]*fx_dot_fx[l]*fx_dot_Dfx[m]*fx_dot_fx[n];
                            re[id]+= pg->pgps[i+l-2][j+m-2][k+n-2].normal_coeff[2]*fx_dot_fx[l]*fx_dot_fx[m]*fx_dot_Dfx[n]; 
                         //   sum++;
                        }
                    } 
                } 
                re[id]/=pg->width; 
            }
        }
    }

    printf("sum:%d\n",sum);
    return re; 
}

